rm(list=ls())
graphics.off()

library(readxl)
library(ggplot2)
library(dplyr)
library(tidyverse)
library(ggpubr)
library(rstatix)
library(nlme)
library(matrixcalc)
library(Biodem)
library(stats)
library(multcomp)
library(hrbrthemes)
library(viridis)
library(gapminder)
library(devtools)
library(stringr)
library(patchwork)
library(gridExtra)
library(cowplot)
library(extrafont)
library(gtable)
library(ggtext)
library(effsize)
library(data.table)
library(psych)
library(psychotools)
library(psychTools)
library(BlandAltmanLeh)
library(pROC)
library(ROCR)

######################### START ###########################################

setwd("//nausers01/User/GOIKOG/Desktop/MatLab/Output/IntensityMeasures3")

data <- read_xlsx("//nausers01/User/GOIKOG/Desktop/MatLab/Output/IntensityMeasures3/MeanRRNormalizedTable.xlsx", col_names = TRUE)


#RRNormalizedTable
RMSSD <- bind_cols(data[,1],data[,24:29])
HR <- bind_cols(data[,1],data[,3:8])


colnames(RMSSD) <- c("ID","aclow","bcchall","cchigh","dmlow","emchall","fmhigh")
colnames(HR) <- c("ID","aclow","bcchall","cchigh","dmlow","emchall","fmhigh")

RMSSD = RMSSD[-c(1,14,18),] #delete bad quality 

total_scores <- RMSSD[,1:7]

RMSSD_long <- pivot_longer(RMSSD, cols=c("aclow","bcchall","cchigh","dmlow","emchall","fmhigh"))
div <- rep(c("Mental","Motor"), each = 3)
RMSSD_long$div <- rep(div,nrow(RMSSD))

RMSSD_long = split(RMSSD_long, RMSSD_long$div)

RMSSD_coglong = RMSSD_long$Mental[,1:3]
  
RMSSD_motlong =  RMSSD_long$Motor[,1:3]

###

plot <- ggplot(RMSSD_long, aes(name, value, ID, fill = name)) +
  geom_boxplot() +
  labs(x = 'Intensity Level',
       title = 'HRV',
       y = "RMSSD (RR Normalized)") +
  theme(plot.title = element_text(hjust = 0.5, size = 20),
        legend.title=element_blank(),
        panel.background = element_blank(),
        panel.grid.major.y = element_line(colour = "grey", size = 0.5),
        axis.title = element_text(size = 13),
        axis.text.x = element_blank(),
        axis.ticks.x = element_blank(),
        strip.text.x = element_text(size = 13)) +
  scale_x_discrete(label = element_blank()) +
  scale_y_continuous(limits = c(-0.015,0.223)) +
  guides(fill="none") +
  scale_fill_manual(values=c("steelblue4", "lightblue1", "cadetblue3", "steelblue4", "lightblue1", "cadetblue3")) +
  annotate(geom = "text", x = 0.83, y = -0.009, label = "Very Low", hjust = 0, vjust = 0, size = 3, color = "gray50") +  
  annotate(geom = "text", x = 0.85, y = -0.015, label = "Intensity", hjust = 0, vjust = 0, size = 3, color = "gray50") +
  annotate(geom = "text", x = 1.85, y = -0.009, label = "Challenging", hjust = 0, vjust = 0, size = 3, color = "gray50") +
  annotate(geom = "text", x = 1.85, y = -0.015, label = "Intensity", hjust = 0, vjust = 0, size = 3, color = "gray50") +
  annotate(geom = "text", x = 2.85, y = -0.009, label = "Very High", hjust = 0, vjust = 0, size = 3, color = "gray50") +
  annotate(geom = "text", x = 2.87, y = -0.015, label = "Intensity", hjust = 0, vjust = 0, size = 3, color = "gray50") +
  facet_grid(. ~ div, scales = "free", space = "free")


#ggsave("HRVplot.png", plot, width = 10.5, height = 6.92)

###Statistics

###Test Normality
Norm.RMSSD.CogEasy = shapiro.test(RMSSD$aclow)
Norm.RMSSD.CogMiddle = shapiro.test(RMSSD$bcchall)
Norm.RMSSD.CogDifficult = shapiro.test(RMSSD$cchigh)

Norm.RMSSD.MotEasy = shapiro.test(RMSSD$dmlow)
Norm.RMSSD.MotMiddle = shapiro.test(RMSSD$emchall)
Norm.RMSSD.MotDifficult = shapiro.test(RMSSD$fmhigh)

###Perform Test and Post-Hoc
### Cognitive 
if (Norm.RMSSD.CogEasy$p.value > 0.05 & Norm.RMSSD.CogMiddle$p.value > 0.05 & Norm.RMSSD.CogDifficult$p.value > 0.05){
  print(anova_test(data=RMSSD_coglong, dv = value, wid=ID, within = name))
  pairwise.t.test(RMSSD_coglong$value, RMSSD_coglong$name, p.adjust.method = 'bonferroni', paired = TRUE)
} else { 
  print(friedman.test(data=RMSSD_coglong, value ~ name|ID))
  pairwise.wilcox.test(RMSSD_coglong$value, RMSSD_coglong$name, p.adjust.method = 'bonferroni', paired = TRUE)
}

###Motor
### Cognitive 
if (Norm.RMSSD.MotEasy$p.value > 0.05 & Norm.RMSSD.MotMiddle$p.value > 0.05 & Norm.RMSSD.MotDifficult$p.value > 0.05){
  print(anova_test(data=RMSSD_motlong, dv = value, wid=ID, within = name))
  pairwise.t.test(RMSSD_motlong$value, RMSSD_motlong$name, p.adjust.method = 'bonferroni', paired = TRUE)
} else { 
  print(friedman.test(data=RMSSD_motlong, value ~ name|ID))
  pairwise.wilcox.test(RMSSD_motlong$value, RMSSD_motlong$name, p.adjust.method = 'bonferroni', paired = TRUE)
}


#ICC analysis

data_TestRR <- read_xlsx("//nausers01/User/GOIKOG/Desktop/MatLab/Output/IntensityMeasures3/TestRRNormalizedTable.xlsx", col_names = TRUE)
data_ReTestRR <- read_xlsx("//nausers01/User/GOIKOG/Desktop/MatLab/Output/IntensityMeasures3/ReTestRRNormalizedTable.xlsx", col_names = TRUE)


## RMSSD RR normalized (used this for the paper)

## Cognitive

### Easy
data_RMSSD_Cog_easy = data.table(data_TestRR$RMSSD_Cog_easy, data_ReTestRR$RMSSD_Cog_easy)
#data_RMSSD_easy_RRn = data_RMSSD_easy_RRn[-3,] #delete unwanted ID (1, 12, 14?)

ICC_RMSSD_Cog_easy =ICC(data_RMSSD_Cog_easy, missing=TRUE, alpha = 0.05, lmer = TRUE, check.keys = FALSE) #last 2 variables not needed, they are so by default

BAP_RMSSD_Cog_easy<- bland.altman.plot(data_TestRR$RMSSD_Cog_easy,data_ReTestRR$RMSSD_Cog_easy, main="RR Normalized Test-Retest RMSSD Easy", xlab="Means", ylab="differences")

print(ICC_RMSSD_Cog_easy)

###Middle
data_RMSSD_Cog_middle = data.table(data_TestRR$RMSSD_Cog_optimal, data_ReTestRR$RMSSD_Cog_optimal)
#data_RMSSD_Cog_middle = data_RMSSD_Cog_middle[-3,] #delete unwanted ID (1, 12, 14?)

ICC_RMSSD_Cog_middle =ICC(data_RMSSD_Cog_middle, missing=TRUE, alpha = 0.05, lmer = TRUE, check.keys = FALSE) #last 2 variables not needed, they are so by default

BAP_RMSSD_Cog_middle <-bland.altman.plot(data_TestRR$RMSSD_Cog_optimal,data_ReTestRR$RMSSD_Cog_optimal, main="RR Normalized Test-Retest RMSSD Optimal", xlab="Means", ylab="differences")

print(ICC_RMSSD_Cog_middle)


###Difficult
data_RMSSD_Cog_difficult = data.table(data_TestRR$RMSSD_Cog_difficult, data_ReTestRR$RMSSD_Cog_difficult)
#data_RMSSD_Cog_difficult = data_RMSSD_Cog_difficult[-3,] #delete unwanted ID (1, 12, 14?)

ICC_RMSSD_Cog_difficult =ICC(data_RMSSD_Cog_difficult, missing=TRUE, alpha = 0.05, lmer = TRUE, check.keys = FALSE) #last 2 variables not needed, they are so by default

BAP_RMSSD_Cog_difficult<- bland.altman.plot(data_TestRR$RMSSD_Cog_difficult,data_ReTestRR$RMSSD_Cog_difficult, main="RR Normalized Test-Retest RMSSD Difficult", xlab="Means", ylab="differences")

print(ICC_RMSSD_Cog_difficult)

## Motor

### Easy
data_RMSSD_Mot_easy = data.table(data_TestRR$RMSSD_Mot_easy, data_ReTestRR$RMSSD_Mot_easy)
#data_RMSSD_easy_RRn = data_RMSSD_easy_RRn[-3,] #delete unwanted ID (1, 12, 14?)

ICC_RMSSD_Mot_easy =ICC(data_RMSSD_Mot_easy, missing=TRUE, alpha = 0.05, lmer = TRUE, check.keys = FALSE) #last 2 variables not needed, they are so by default

BAP_RMSSD_Mot_easy<- bland.altman.plot(data_TestRR$RMSSD_Mot_easy,data_ReTestRR$RMSSD_Mot_easy, main="RR Normalized Test-Retest RMSSD Easy", xlab="Means", ylab="differences")

print(ICC_RMSSD_Mot_easy)

###Middle
data_RMSSD_Mot_middle = data.table(data_TestRR$RMSSD_Mot_optimal, data_ReTestRR$RMSSD_Mot_optimal)
#data_RMSSD_Mot_middle = data_RMSSD_Mot_middle[-3,] #delete unwanted ID (1, 12, 14?)

ICC_RMSSD_Mot_middle =ICC(data_RMSSD_Mot_middle, missing=TRUE, alpha = 0.05, lmer = TRUE, check.keys = FALSE) #last 2 variables not needed, they are so by default

BAP_RMSSD_Mot_middle <-bland.altman.plot(data_TestRR$RMSSD_Mot_optimal,data_ReTestRR$RMSSD_Mot_optimal, main="RR Normalized Test-Retest RMSSD Optimal", xlab="Means", ylab="differences")

print(ICC_RMSSD_Mot_middle)


###Difficult
data_RMSSD_Mot_difficult = data.table(data_TestRR$RMSSD_Mot_difficult, data_ReTestRR$RMSSD_Mot_difficult)
#data_RMSSD_Mot_difficult = data_RMSSD_Mot_difficult[-3,] #delete unwanted ID (1, 12, 14?)

ICC_RMSSD_Mot_difficult =ICC(data_RMSSD_Mot_difficult, missing=TRUE, alpha = 0.05, lmer = TRUE, check.keys = FALSE) #last 2 variables not needed, they are so by default

BAP_RMSSD_Mot_difficult<- bland.altman.plot(data_TestRR$RMSSD_Mot_difficult,data_ReTestRR$RMSSD_Mot_difficult, main="RR Normalized Test-Retest RMSSD Difficult", xlab="Means", ylab="differences")

print(ICC_RMSSD_Mot_difficult)



